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ABSTRACT 

We study the type III migration of a Saturn mass planet in low viscosity discs. The 
planet is found to experience cyclic episodes of rapid decay in orbital radius, each 
amounting to a few Hill radii. We find this to be due to the scattering of large- 
scale vortices present in the disc. The origin and role of vortices in the context of 
type III migration is explored. It is shown through numerical simulations and semi- 
analytical modelling that spiral shocks induced by a sufficiently massive planet will 
extend close to the planet's orbital radius as well as being global prominent features. 
The production of vortensity across shock tips results in thin high vortensity rings with 
a characteristic width of the local scale height. For planets with masses equal to and 
above that of Saturn, the rings are co-orbital features extending all the way around the 
orbit. Linear stability analysis show such vortensity rings are dynamically unstable. 
There exists unstable modes that are localised about local vortensity minima which 
coincide with gap edges. Simulations show that vortices are non-linear a outcome. 

We used hydrodynamic simulations to examine vortex-planet interactions. Their 
effect is present in discs with kinematic viscosity less than about an order of magnitude 
smaller than the typically adopted value of v — 10~ 5 £7 p r p (0) 2 , where r p (0) and ft p 
are the initial orbital radius and angular velocity of the planet respectively. We find 
that the magnitude of viscosity affects the nature of type III migration but not the 
extent of the orbital decay. The role of vortices as a function of initial disc mass is 
also explored and it is found that the amount of orbital decay during one episode 
of vortex-planet interaction is independent of initial disc mass. We incorporate the 
concept of the co-orbital mass deficit in the analysis of our results and link it to the 
presence of vortices at gap edges. 
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1 INTRODUCTION 

The importance of planet migration was realised with the 
disco very of the first exopla net with an orbital period of 4 
days IjMavor fc Queloal 19951 ). Such planets, classified as 'hot 
Jupiters' orbit so close to their host stars it is difficult to 
understand their formation in situ. It is thought they form 
further out then migrate inwards due to torques from the 
gaseous protoplanetary disc. 

Tidal interactions between a protoplanet and proto- 
planetary disc were in fact studied well befor e obse rva- 
tions of exoplanets (|Goldreich fc Tremaine]|l979l . [198(31 ). In 
these early studies, the protoplanet is treated as a small 
perturbation, the linearised hydrodynamic equations are 
then solved for the disc response. Density waves launched 
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at Lindblad resonances generally lead to decay of the 
planet's semi-ma jor axis. T his scenario is referred to as 
type I migrati on llWa rd 1997). Rec ent semi-analytical treat- 
ments include Tanaka et al.l (i2002) for isothermal discs and 
iPaardekooper et alj 1 20091 ) for a non-isothermal equation of 
state. For planet masses above that of Jupiter, linear theory 
ceases to apply and the non-linear disc response produces 
a sur face density gap about the planet (|Lin fc Papaloizovj 
1986). In this migration mode, now called type II, the 
planet's orbital radius is locked with the disc viscous evolu- 
tion. That is, it drifts towards the central star along with 
disc material, while residing inside the gap. 

Recently, a new migration mod e, called runaway or 
type I II migration, was introduce d bvlMasset fc Papaloizoul 
(l2003h and furthe r discussed by lArtvmowicz ( 2004a ) and 
iPapaloizoul (|2005l ). Type III migration is a self-sustaining 
mechanism with short migration timescales O(10 2 ) or- 
bits. A series of detailed numerical studies was presented by 
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iPepliriski et~al1 (|2008al lbll3) , who considered a Jupiter-mass 
planet and included corrections for self-gravity and thermal 
effects. In their standard case, the planet's semi-major axis 
is reduced by a factor of ~ 3 within 50 orbital periods. They 
also concluded that type 111 migration is strongly dependent 
on the flow inside the Hill radius of the planet. 

Type III migration is also applicable to intermediate 
mass planets comparable to Saturn in massive discs. The 
planet opens a partial gap which, in the case of a migrating 
planet, allows a net drift of disc material across the planet's 
orbital radius through executing U-turns at the end of horse- 
shoe orbits. The change in the fluid elements' orbital radius 
implies a torque on the planet that is proportional to its 
migration rate. The migration rate also increases with the 
co-orbital mass deficit Sm, which is proportional to the dif- 
ference between surface density of the co-orbital region and 
that of the orbit-crossing flow. 

Most previous disc-planet simulations have included 
large enough viscosity to suppress instabilities resulting 
in a smooth behaviour. Even numerical viscosity is suf- 
fic ient to suppres s instab ilities. For example, the studies 
of IPepliriski et all f2008a,b,c) did not display instabilities 
even with zero physical viscosity. Howev er, they employed a 
Carte sian code which, as pointed out bv lde Val-Borro et al.l 
(120071) . is probably too diffusive to allow steep gradients to 
form and become unstable. 

In this work, we study a modified form of type III mi- 
gration, which is non-smooth, induced by large-scale vortices 
originating from instabilities. The vortex forms outside the 
Hill sphere and flows across the co-orbital region. Chang- 
ing the planet softening length or adopting the equation of 
state of IPepliriski et all (|2008al ). which modifies flow inside 
the Hill sphere, is found not to significantly affect the vortex- 
planet interaction described in this paper. 

The linear theory of instabilities in inviscid accre- 
tion discs associated with ring str uctures, sharp edges 



2 BASIC EQUATIONS AND MODEL 

We consider a planet of mass M p orbiting a central star of 
mass M*. We adopt a cylindrical coordinate system (r, (f>, z) 
where z is the vertical coordinate increasing in the direction 
normal to the disc plane for which the unit vector is k. 
We integrate vertically to obtain a two-dimensional flat disc 
model, for which the governing hydrodynamic equations in 
a frame uniformly rotating with angular velocity fi p k are 
the continuity equation 

§ = -£V.u, (i) 

and the equation of motion 

^ + 2n„kAu = -ivp-v$ cff , (2) 

where D/Dt is the total derivative and the vertically inte- 
grated pressure P = c 2 (r)E with c 3 = /i(GM*/r) 1//2 . Here 
E is the surface density, u is the velocity and the effective 
gravitational and centrifugal potential is given by 
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(Massct 2002). The last two terms on RHS are indirect terms 
accounting for acceleration of the primary. In the above, the 
cylindrical coordinates of the planet that is confined to re- 
main in the disc plane are (r p , <f> p , 0) and the softening length 
e = 0.6H(r p ) where H (r) = hr is the disc semi-thickness and 
h — 0.05 the constant aspect ratio. All numerical and an- 
alytical work are based on the two-dimensional equations, 
but H and e accounts for vertical structure. 
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20011 ). In the context of protoplanetary discs, steep gradi- 



ents aris e from gap edges associate d with sufficiently massive 
planets. Ide Val-Borro et"al] (I2007T ) performed simulations of 
a Jupiter-mass planet held on fixed orbit to show the for- 
mation of large-scale vortices at gap e dges. The migrating 
case was considered bv lOu et al.l ((2007) for a Neptune-mass 
planet. The authors find non-smooth migration associated 
with non-axisymmetric surface density enhancements near 
the gap edge. In this paper we examine this effect in more 
detail for a Saturn mass planet undergoing type III migra- 
tion. 

This paper is organised as follows. Section [2] describes 
the set of equations that we used to model the disc-planet 
system and in Section[3]we describe the formation of vorten- 
sity rings using numerical simulations and a semi-analytical 
model. We then study their dynamical stability in Section 
|3] In Section \5\ we present hydrodynamic simulations of type 
III migration as a function of viscosity, high-lighting the ef- 
fect of vortices at low viscosity. We analyse the inviscid case 
in detail in Section refinviscid and consider varying the disc 
mass and planet mass. Finally in Section[7]we summarise 
our results. 



The vortensity being the ratio of the z component of the 
vorticity (hereafter just called the vorticity, ui) to the surface 
density is defined to be 

V = ^ = ^ — -, (4) 



UJ 

E = 
where to. 



z ■ V A u is the relative vorticity seen in the ro- 
tating frame. It is well known that In barotropic flows with- 
out shocks, it follows from ([TJ and © that the vortensity 
rj = uj/T, is conserved for a fluid particle. When an isother- 
mal equation of state with variable sound speed as is used 
here is adopted, vortensity is no longer strictly conserved. 
However, the sound speed varies on a global scale so that 
when phenomena are considered on a local scale, vorten- 
sity is conserved to a good approximation in the absence of 
shocks. 



2.2 Hydrodynamic simulations 

Our work is based on numerical simulations and specific 
parameters depend on the problem at hand, as is described 
in the following sections. Here we present the general set up. 

For convenience we adopt dimensionless units such that 
the orbital radius of the planet, r p (t), is initially r p (0) = 
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2. Time is expressed in units of the initial orbital period 
Po = 2n/Q, p , where Q p = ^/GM»/r|(0) is the planet's initial 
Keplerian angular velocity. The unit of mass is taken to be 
the central mass M, which for working purposes may be 
considered to be IMq. The disc has an initially uniform 
surface density E = Eo x 10~ 4 . Some simulations include 
a constant kinematic viscosity v = Vq X 10 -5 in physical 
units of r5(0)f2 p , with vo being a dimensionless constant. 
We considered Eo = 1 — 9 and vq = — 1. 

We use th e FARGO code to solve for the disc response 
|Masse t 2000a. b). FARGO evolves the hydrodynamic equa- 
tions in global cylindrical co-ordinates centred on the central 
object. We used a non-rotating frame for the simulations 
described here. Indirect terms arising from the non-inertial 
frame are included in the potential. The code uses a finite- 
difference scheme with van Leer upwind a dvection , similar 
to the ZEUS code (|Stone fc Normanll 19921 ) but it employs a 
modified algorithm for azimuthal transport that allows for 
large time steps. We take 27r-periodic boundary condition 
in azimuth and wave damping boundary conditions at disc 
boundaries (|de Val-Borrdl2006h . We remark that vortensity 
ring formation occurs as a result of shocks near the planet 
and therefore is unaffected these boundary conditions. 

In the case where the planet is allowed to migrate, a 
fifth order Rungke-Kutta method was used to integrate its 
equation of motion. Simulations ran with halved time- 
step show very similar results to the corresponding 
standard set-up. In particular, the extent of planet 
migration by vortex-scattering and the number of 
events of scattering is the same. We believe the 
RK5 integrator is sufficiently accurate to capture 
this interaction, which can be regarded as a two- 
body problem. 



3 RING STRUCTURES IN DISC-PLANET 
INTERACTIONS 

In this section we present hydrodynamic simulations to show 
that narrow surface density rings are brought about as a 
consequence of the fact that highly peaked vortensity (be- 
ing the ratio of the vorticity to the surface density) rings are 
produced by flow through quasi-steady shocks located close 
to the planet. For sufficiently massive planets, the associ- 
ated vortensity generation occurs as fluid elements execute 
a horseshoe turn in the co-orbital region. Focusing on such 
cases, we construct a simple model that enables an esti- 
mate of the shock location to be made together with rate 
of the vortensity generation as material flows through it. 
Our model is in good agreement with numerical simulations 
and confirms the process of vortensity generation within a 
planet's co-orbital region which later plays an important role 
in modifying the flow structure there. 



3.1 Vortensity generation by the shocks induced 
by a Saturn mass planet 

We first consider the case with M p = 2.8 x 10~ 4 M, cor- 
responding to a Saturn mass planet around a solar mass 
star. Here, the planet potential is switched on over 5Po- The 
computation domain is r = [1,3] and the resolution was 



N r x N$ — 1024 x 3072, corresponding to radial and az- 
imuthal grid spacings of Ar ~ 0.02rh and rA<fi ~ 0.05rh re- 
spectively, where = (M p /(3M*)) 1 ^ 3 r p is the planet's Hill 
radius. The planet is held on fixed circular orbit and the disc 
has density Eo = 1 with no explicit viscosity (v — 0). 

Figfl] shows the vortensity field at t — 7.07 close to 
the planet. Vortensity is generated/destroyed as material 
passes through the two spiral shocks. For the outer shock, 
vortensity generation occurs for fluid elements executing a 
horseshoe turn (\r — r p \ < Th) while vortensity is reduced for 
fluid that passes by the planet, but the change is smaller in 
magnitude in the latter case. The situation is similar for the 
inner shock, but some post-shock material with increased 
vortensity continues to pass by the planet. This stream be- 
gins at r ~ — Th but such feature is absent at the outer 
shock, suggesting a lack of symmetry about r — r p , possibly 
resulting from the non-uniform vortensity background be- 
ing oc r -3 / 2 . Note that although a pre-shock fluid element 
that would be on a horseshoe trajectory, may in fact pass by 
the planet after crossing the shock, it is clear that vortensity 
rings originate from passage through shock fronts interior to 
the co-orbital region that would correspond to the horseshoe 
region for free particle motion. 

The streams of high vortensity eventually move around 
the whole orbit outlining the entire co-orbital region. Fig. 
[1] shows that they are generated along a small part of the 
shock front of length ~ r^. This results in thin rings with 
a similar radial width. The fact that they originate from 
horseshoe material can enhance the contrast as post-shock 
inner disc horseshoe material is mapped from r — x to r+x to 
become adjacent to post-shock outer disc material passing 
by the planet. 

We also show in Fig. f2] long-term evolution of average 
co-orbital properties from a corresponding low resolution 
run (N r x A^ = 256 x 768). The co-orbital region is taken 
to be the annulus [r p — x s ,r p + x s ]. We fix x s — 2.5rh, 
as is typically mea s ured from hydrodynamic simulat ions 
i|Artvmowiczl l2004bl : iPaardekooper fc Papaloizoul I2009T ) for 
intermediate or massive planets. In Appendix |X] we show 
that in the particle dynamics limit x s < 2.3rh, comparable 
to the value adopted above. 

Vorticity generation occurs within t < 25 orbits, after 
which it remains approximately steady. For a Jupiter-mass 
planet the time taken to reach this state is about 50 or- 
bits, but most of the vorticity generation also occur in the 
first ~ 25 orbits. It is important to note that subsequent 
vortensity increases is in narrow rings and fluctuations are 
due to instabilities associated with the rings. The average 
surface density falls more gently as the planet opens a gap, 
which is a requirement for the rings to be self-supported. 
Consequently, we observe a rise in co-orbital vortensity. Fig. 
[2] reflects modification of co-orbital properties on dynamical 
time-scales due to shocks, so we expect migration mecha- 
nisms which depend on co-orbital structure to be affected. 

3.2 Location of the vortensity generation region 
for different planet masses 

The process of formation of vo rtensity ring s disc ussed 
here diff e rs fro m that observed in iKoller et alj (|2003T l and 
iLi et all (120051 ), where the rings are generated by shocks 
outside the co-orbital region. This is because the authors 
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7.07 orbits, r p =2.00 
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Figure 1. Vortensity generation and destruction across shocks 
induced by a Saturn-mass planet in an inviscid disc. The plot 
shows lnr; = ln(o;/E). Regions of increased vortensity are clearly 
visible as half horseshoes above and below the planet. The in- 
crease occurs as material passes through parts of the shock fronts 
that extend into the co-orbital region. 
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Figure 2. Long-term evolution of co-orbital properties for a 
Saturn-mass planet held on fixed orbit. Vorticity, density and 
vortensity perturbations are expressed relative to their values at 
t = 0. Angle brackets denote an average over the co- 
orbital region defined by the annulus r = [r p — 2.5ri l ,r p + 
2.5x s ]. 

used a smaller planet mass. To illustrate the effect of reduc- 
ing the mass, we ran simulations with N r x = 256 x 768 
for M p x 10 4 = 1, 2.8 and 10. Fig. [3] compares azimuthally 
averaged vortensity perturbations induced in each case after 
14.14 orbits. The intermediate and high mass cases are qual- 
itatively similar, having A(o;/E) maxima at r — r p ~ ±2r; t 
and minima at r — r p ~ 3r^. The magnitude of A(cj/E) in- 
creases with M p because higher masses induce higher Mach 
number shocks. As the half-width of the horseshoe region is 
x s = 2.5rfc for such masses, vortensity rings are co-orbital 
features. 

The low mass M p = 10~ 4 case has much smaller 
|A(cj/E)|. There is a vortensity maximum at r — r p = 
— 3r^ but nothing of similar magn itude at r — r p > 
0. IPaardekooper &: Papaloizoul (|2009l ) found the co-orbital 
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Figure 3. Azimuthally averaged vortensity perturbation for dif- 
ferent planet masses. 

half-width, in the limit of zero softening for low mass plan- 
ets, to be: 

x a = im{M p /h) 1/2 r p . (5) 

For M p = 10~ 4 , h = 0.05, x s = 2.33r h . Non-zero soften- 
ing gives a smaller x s - Hence, vortensity rings for low mass 
planets occur outside the co-orbital region as is confirmed 
in our simulations and they are very much weaker. 

3.3 A semi-analytic model for co-orbital 
vortensity generation by shocks 

We now construct a semi-analytic model describing the pro- 
cess of shock-generation of vortensity. In this way we gain an 
understanding of the physical processes involved and later 
how they lead to strong torques and fast migration of the 
planet and under what conditions simulations can represent 
them accurately. 

More specifically we model the outer spiral shock in Fig. 
[1] To do this we need the pre-shock flow field, the shock 
front location and then to evaluate the vortensity change 
undergone by material as it passes through the shock. We 
now consider these in detail. 



3.3.1 Flow field 

We first simplify the geometry by adopting the shearin g box 
approximation (e.g. IPaardekooper fc Papaloizoul I2009T ). As 
we are concerned with the flow near the planet, we consider 
a local Cartesian co-ordinate system (x,y) co-rotating with 
the planet at angular speed Q p and with origin at its centre 
of mass. Without the planet, the velocity field is Keplerian 
u = (0, — 3Q p x/2). In order to deal with the pre-shock flow 
we make the assumption that pressure effects can be ig- 
nored when compared to the planet potential. This ballistic 
approximation is appropriate for a slowly varying supersonic 
flow as is expected to be appropriate for the pre-shock fluid. 
The velocity u satisfies the local form of equation ([2]) with 
P — and the effective potential 

GMp 3 2 2 / c n 
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and thus follows from particle dynamics. The indirect terms 
are neglected in this approximation. We write a fluid parti- 
cle's trajectory in a steady state flow in the form of x — x(y) 
and the corresponding velocity field as u = u(y). Noting 
that on a particle trajectory we have -jjj = u v£j, it follows 
from the local form of Q that on a particle trajectory we 
have the following system of three simultaneous first order 
differential equations: 



du 



dy 
d 



2M p y 



(x 2 + y 2 + e 2 ) 3 / 2 



(UxUy) = Uy 



M p x 



(a; 2 



_|_ £ 2)3/2 



+ U X Q 



(xuy) = xQ + u x u y . 



(7) 



(8) 
(9) 



dy 



dy 

We use these to solve for the state vector \J(y) = 
[u^, u x v 2 , xuy] in x > for a particular particle. The bound- 
ary condition is 

(u y , u x , x) — > (— 3Q p xo/2, 0, xo) as y 



(10) 

where x — xo is the particle's unperturbed path. The totality 
of paths obtained by considering a continuous range of xo 
then constitutes the flow field. Having obtained the velocity 
field, we use vortensity conservation to obtain E. The surface 
density is then 



Z(x,y) = ^(u r + 2n p ) 



(11) 



where the unperturbed absolute vorticity is u = f2 p /2. Nu- 
merically, we calculated the relative vorticity by relating it 
to the circulation through 

Ur ~ ~Ka / u ' dl (12) 

where the integration is taken over a closed loop about the 
point of interest and AA is the enclosed area. This avoids 
errors due to numerical differentiation on the uneven grid in 
(x, y) generated by solving the above system. 

3.3.2 Location of the shock front 

The physical reason for shock formation is the fact that the 
planet presents an obstacle to the flow. Where the relative 
flow is subsonic, the presence of this obstacle can be commu- 
nicated to the fluid via sound waves emitted by the planet. 
In regions where the relative flow is supersonic, the fluid 
is unaware of the planet via sound waves (but the planet's 
gravity is felt). We estimate the location of the boundary 
separating these two regions by specifying an appropriate 
characteristic curve or ray defining a sound wave. This is 
a natural location for shocks. App lying this idea to Keple- 
rian flow, IPapaloizou et all ([2004 ) obtained a good match 
between the predicted theoretical shock front and the wakes 
associated with a low mass planet. For general velocity field 
u, the characteristic curves satisfy the equation 



dy s 
dx 



u x u y — \Ju\+ii\ 



(13) 



where ft; = Ui/c B . The positive sign of the square root 
has been chosen so that the curves have negative slope in 



the domain x > with u y < 0. As the fluid flows from 
super-sonic (y > y s ) to sub-sonic (y < y s ), fluid located at 
y = y s begins to know about the planet through pressure 
waves (|Papaloizou et alj|2004 ). In Keplerian flow, the rays 
defining the shock fronts originate from the sonic points at 
x — ±2H/3, y — 0. In a general flow, sonic points where 
(|u| = c s ), at which the rays may start, lie on curves and 
can occur for x < 2hr p /3. The starting sonic point for solv- 
ing equation ( I13p that we eventually adopted has x = and 
the value of y > that is furthest possible from the planet. 



3.3.3 Vorticity and vortensity changes across a shock 

The jump in absolute vorticity [u] is readily obtained by 
resolving the flui d motion paralle l and perpendicular to the 
shock front (eg. iKevlahar] 1 19971 ) . As we do not solve the 
energy equation, shocks are l ocally isotherma l and our ex- 
pression differs from those of iKevlahanl l| 19971 ). accordingly 
a brief derivation of [lj] is presented in Appendix B. The 
result for a steady shock is 



(M 2 - l) 2 du ± 
dS~ 



M 2 



+ (M - 1)oj- 



M 2 



dS' 



(14) 



where u± is the pre-shock velocity component perpendicu- 
lar to the shock front, M = u±/c a is perpendicular Mach 
number and c a = hr~ x ^ 2 is the sound speed. d/dS is the 
derivative along the shock. It is important to note that for 
()14|) to hold, the direction of increasing S, the direction of 
positive itx, and the vertical direction should form a right 
handed triad. We take increasing S as moving away from 
the planet. The local isothermal equation of state produces 
the last term on RHS of equation (|14|l . Because of the slow 
variation of c s its contribution is not important in this ap- 
plication. 

The vortensity jump [w/E] follows immediately from 
equation (|14[) as 



ffi] 



(M 2 - l) 2 du A 



EM 4 dS \T,M 2 u 



f M 2 -1\ 9c? 

dS 1 



(15) 



which e s sentia lly reduces to the expression derived by 
iLi et al.l (120051 ) if c s = constant (and a sign change due 
to different convention). The sign of [oj/E] depends mainly 
on the gradient of u± (or M) along the shock. As in our case 
mi < 0, the width of the increased vortensity rings produced 
is determined by the length along the shock where \M\ is in- 
creasing. Note that [w/E] does not depend explicitly on the 
pre-shock vortensity, unlike the absolute vorticity jump. 



3.3.4 Comparison of the results of the semi-analytic 
model with numerical simulations 

We now compare the results of hydrodynamic simulations 
with those obtained from the model. First we illustrate the 
particle paths that constitute the flow field together with 
the shock fronts obtained by assuming coincidence with 
the characteristic curves that are obtained from the semi- 
analytic model in Fig. [4] A polynomial fit to the simu- 
lation shock front is also shown. Particle paths cross for 
x > rh, y < —2rh so that the neglect of pressure becomes 
invalid. Accordingly the solution to equation Q should not 
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7.07 orbits, M =2.8x10" 




4.0 



Figure 4. Solid lines: particle paths from the reduced zero- 
pressure momentum equations (equations [7] - [5J); thick lines: 
curves composed of sonic points on which [u| = c s ; dotted lines: 
theoretical shock fronts; dash-dot line: solution to equation J 13D 
for Keplerian flow; dashed line: polynomial fit to simulation shock 
front. The actual shock front begins at a sonic point around 
x = 0.2rv 



be trusted in this region. Fortunately, vortensity generation 
occurs at a distance from the planet that is within rn, where 
pre-shock particle paths do not cross. 

In Fig. [J the estimated shock location is qualitatively 
good and tends to the Keplerian solution further away. The 
important feature is that the shock can extend close to 
x = 0, across horseshoe orbits. If the flow were purely Keple- 
rian there could be no significant vortensity generation close 
to x — because the flow becomes sub-sonic there. Further- 
more, only circulating fluid can be shocked in that case. 
Shock-generation of vortensity inside the co-orbital region 
is only possible for sufficiently massive planets that induce 
large enough non-Keplerian velocities. 

The key quantity determining vortensity generation is 
the perpendicular Mach number. Fig. [5] compares M from 
simulation and model. Although the semi-analytic model 
gives a shock Mach number that is somewhat smaller than 
found from the simulation, all curves have \M\ increasing 
from x — — > 1.37-^ which is the important domain for 
vortensity generation. Thus equation (I15[) implies vorten- 



sity generation in all cases. Fig. 6(a) illustrates various com- 
binations of semi-analytic modelling and simulation data 
that have been used to estimate the vortensity jump across 
the shock. The qualitative similarities between the various 
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Figure 5. Perpendicular Mach number squared M 2 along the 
outer spiral shock illustrated in Fig.[T] Asterisks: simulation data; 
dotted line: M 2 obtained from the simulation shock and the 
semi-analytic flow; solid line: M 2 obtained from the semi-analytic 
shock front and flow. 



curves confirms that vortensity generation occurs within co- 
orbital material about one Hill radius away from the planet. 
It is shocked as it undergoes horseshoe turns. 

Assuming that material is mapped from x — > — x as it 
switches to the inner leg of its horseshoe orbit, we expect 
the outer spiral shock to produce a vortensity ring peaked 
at r — r p ~ — 0.5rh of width 0(rh) (= O(H)), with a similar 
discussion applying to the inner shock which is expected 
to produce a vortensity ring peaked at r — r p ~ 0.5r"h . 
Of course, as a fluid element moves away from the U-turn 
region, r — r p \ increases, but it remains on a horseshoe orbit. 
Thus, thin vortensity rings are natural features of the co- 
orbital region for such planet masses. 

We have also tested our model for M v 



Fig. 6(b) There is still a good qualitative match between 



the simulation and the model; even though lowering M v 
makes the zero-pressure approximation, adopted to deter- 
mine the semi-analytic flow field, less good. Decreasing M v 
shifts vortensity generation away from the planet in the 
semi-analytic model, as is also observed in the hydrodynamic 
simulation (Fig. [3}. In this case, there is no vortensity gener- 
ation in r — r p < 0.5rh but vortensity rings are still co-orbital 
(with peaks at ~ 0.7rjj). 

We remark that steep vortensit y gradients are associ- 
ated with dynamical instabilities (e.g. lPapaloizou fc Pringlel 
Il985l ) and this is explored in more detail below. 



4 DYNAMICAL STABILITY OF VORTENSITY 
RINGS 

Having established the origin of the vortensity rings and the 
mechanism for producing them, we go on to study the linear 
stability of the shock-modified protoplanetary disc model 
described above. This is an important issue as instability 
can lead to their breaking up into mobile non-axisymmetric 
structures which can affect the migration of the planet sig- 
nificantly. 

The linear stability of inviscid barotropic axisymmct- 
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Figure 6. Semi-analytic and actual vortensity jumps across a 
shock in the co-orbital region. Asterisks are measured from the 
simulations. The dashed lines were obtained by using pre-shock 
simulation data coupled with the jump condition specified by 
cquation( I15| l. The dotted lines were obtained using the semi- 
analytic flow field together with the location of the shock front 
obtained from the simulation. The solid curves correspond to the 
semi-analytic model for both the flow field and shock front. 



ric rings and discs without a planet is well develope d 
jPapaloizou fc Pringld[l984l , Il985l ; iPapaloizou fc LirJI 19891 ) . 
The extension to a non- barotropic equation of state was un - 
dertaken more recently (|Lovelace et al.|[l999l ; lLi et alfe OQO). 
This work indicates that sharply peaked surface density or 
vortensity features in the disc are associated with dynamical 
instabilities. It is directly relevant to the types of configura- 
tion that we have found to be produced b y the disc planet in- 
teract ions described above. For example Ide Val-Borro et al.l 
(120071 ) considered the linear instability of gap edges associ- 
ated with a Jupiter mass planet and connected it to vortex 
formation there. 

Here we consider the stability of the partial gap opened 
by a Saturn- mass planet, for which type III migration is ex- 
pected in a massive disc. In order to be able to do this we 
need to be able to define an appropriate background equilib- 




Figure 7. A contour plot of ln(S/o;) when the vortensity ring 
like structures have formed in the co-orbital region. These are 
almost axisymmetric apart from in the region very close to the 
planet. The centre of the small square marks the planet position. 



rium axisymmetric structure to perturb. We show that this 
is possible. 



4.1 Basic background state 

The basic state should be axisymmetric and time indepen- 
dent with no radial velocity (d/dcj) = d/dt — u r — 0). To 
set this up we suppose the vortensity profile r](r) is known 
(e.g. via shock modelling as above) and that accordingly we 
have 



The radial momentum equation gives 



ldP 
E~dr 



+ 



GM, 



(16) 



(17) 



Using these together with the locally isothermal equation of 
state, P = h 2 GM*Y,/r, h being the constant aspect ratio, 
we obtain a single equation for E in the form 



! lnE 



dr 2 



+ 



2E7? 

vum: 



1-h 2 



+ h 2 



,dlnE 



dr 



1/2 



2h 2 dlnE 
r dr 



(18) 



We solve equation (|18|l for E with rj{r) taken as an azimuthal 
average from the fiducial simulation described in f £|3.1[) at a 
time at which vortensity rings have developed. These struc- 
tures are essentially axisymmetric apart from in the close 
neighbourhood of the planet. They are illustrated in Fig. [7] 
The boundary conditions is E = Eo = 1 at r = 1.1, and 
r = 3. These conditions are consistent with the fact that 
shock-modification of the surface density profile only occurs 
near the planet (r = 2). 

A comparison between hydrostatic surface density and 
angular velocity profiles obtained by solving equation (|18[1 
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with those obtained as azimuthal averages from the cor- 



' 9.45 orbits 



responding simulation is illustrated in Fig. 8(a) and Fig 



8(b) The agreement is generally very good indicating that 
the adoption of the basic axisymmetric state defined by the 
simulation vortensity profile, together with equation (|17[) . 



for stability analysis should be a valid procedure. Fig. 8(a) 
shows that vortensity rings reside in the horseshoe region 
just inside the gap. At a surface density extrema equation 
(|18[) implies that 



Ir 



dr 2 



+ 2Er/ 



1 - h 2 



GM,r 

Since h <C 1, if 77 is sufficiently large then d 2 E/dr 2 > 0. 
Hence vortensity maxima will coincide with surface density 
minima. In our basic state, local minima and maxima are 
separated by O(H) (~ 0.1). As vortensity rings originate 
from spiral shocks, equation ( I18|l demonstrates the link be- 
tween shocks and gap formation. Given the double-ringed 
vortensity profile r/(r), which may be estimated by modelling 
shocks as described above, we can solve equation (|18|) for the 
axisymmetric surface density profile E(r). We will then find 
that in order for the rings to be in hydrostatic equilibrium, 
a gap in the surface density must be present around the 
planet's orbital radius, which is in between the vortensity 
peaks. Sufficiently massive planets induce strong shocks, im- 
plying larger vortensity maxima, and hence deeper surface 
density minima or gaps. 
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4.2 Linearised equations 

Having established the ring-basic state, we perform linear 
analysis to determine its stability. The analysis is for a 
two-dimensional disc and isothermal perturbations. This is 
much simpler than the non-barotropic 2D disc studied by 
iLi et al.l (|2000l ) . who considered adiabatic perturbations and 
hence involved the energy equation, omitted in the present 
work. This is so that linear analysis is consistent with hy- 
drodynamic simulations prese nted above and later. The 3D 
barotropic disc was studied by IPapaloizou fc Pringlei (|l987l ) 
where edges and the role of the vortensity were first consid- 
ered. We obtain the governing equation for locally isother- 
mal perturbations by linearising the continuity equation and 
the equation of motion as seen in the non-rotating frame 

^-+V-(Eu) = (19) 

^ + u-Vu = -ivP-V$, (20) 

where P = c 2 E with c 2 = h 2 GM*/r. The total potential $ 
is assumed fixed. We set 

/E\ /E\ /*E(r)\ 

I u r I — > + Su r (r) x exp i(at + m(f>) 

\U4,J \U4,J \5 Ult> (r)J 

where, the first term on the right hand side corresponds to 
the basic background state, a is a complex frequency and 
m, the azimuthal mode number, is a positive integer. The 
linearised equation of motion (equation (120[) 1 gives 

„2 



K z — G 
2 



._dW 2imQ,W \ 
dr r J 



dW ma 
E?7— 1 W 

dr r 



(21) 



(22) 



(b) n/a k 

Figure 8. A comparison of the axisymmetric surface density and 
angular velocity (scaled by Keplerian speed) profiles obtained by 
solving equation {18}, with the same quantities obtained by az- 
imuthally averaging the simulation data over 2ir. The azimuthally 
averaged vortensity profile from simulation data is also shown 
(dotted). 

The co-orbital region is r = [1.77, 2.22], 

where W = <5E/E is the fractional density perturbation, 
fit 2 = 2fffi77 is the epicycle frequency expressed in terms of 
the vortensity rj, and a = a + mQ(r) = or + mf2(r) + 17 is 
the Doppler-shifted frequency with or and 7 being real. 

Substituting equation (|2ip — \12\ into the linearised 
continuity equation 

1 d 



iaW 



(rT,8u r 



rE dr r 
yields a governing equation for W of the form 

d_( E dW s 
dr 



(23) 



+ 



K — 

m d 
o dr 



dr 



rrj(n 2 



rE 



m 2 E 



GM-tb? 



W = 0. 

(24) 

This equation leads to an eigenvalue problem for the com- 
plex eigenvalue a. Co-rotation resonance occurs at a = 
which requires 7 = for it to be on the real axis. Then for 
the equation to remain regular, the gradient of the terms 
in square brackets must vanish at co-rotation. This results 
in the requirement that the gradient of rjr vanish there. Be- 
cause the sound speed varies with radius, this is slightly dif- 
ferent from the condition that the gradient of rj should van- 
ish which applies in the barotropic strictly isothermal case 
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jPapaloizou fc Pringld[l984l , 1 19851 ; iPapaloizou fe Linll 19891) . 
However, because r\ varies vary rapidly in the region of inter- 
est and the modes we consider are locally confined in radius, 
this difference is of no essential consequence. Lindblad res- 
onances occur when k 2 — o 2 = 0, but as is well known, and 
can be seen from formulating a governing equation for 5u r 
rather than W, these do not result in a singularity. 

4-2.1 Simplification of the governing ODE 

It is useful to simplify equation f|24[) to gain further in- 
sight. To to this we consider modes localised around the co- 
rotation circle such that the condition k 2 3> \o 2 \ is satisfied. 
Beyond this region the mode amplitude is presumed to be 
exponentially small. Now, recognising that GM*h 2 jr = c 2 , 
the ratio of the second to last to last term in equation Q24p 
is 



m 2 c 2 m 2 h 2 

For a thin disc h -C 1, so if we consider m = O(l) this ratio 
is large and the last term in equation (I24|l can be neglected. 
This is also motivated by the fact that only low m modes are 
observed in simulations. Doing this and replacing k 2 — \a 2 \ 
by ft 2 , equation (|24[) reduces to the simplified form 

rf^)^7^1-^K = o. (25) 

dr \ k 2 dr J [ a dr |_ r/ J J 

We comment that in this form equation [25] is valid for any 
fixed c s profile. 
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(b) Outer edge 

Figure 9. m = 3 eigenmodes obtained with the boundary condi- 
tion W = 0. The perturbed quantities are plotted as \W\ (solid, 
black), |<5u r | (dotted, blue), \Su ( j > \ (dashed, green) and scaled by 
their maximum values in r = [1.1,3.0]. The background vorten- 
sity profile is also shown (dash-dot, red), ro = 1.7, 2.3 correspond 
to the inner and outer vortensity minima, respectively. 



4-2.2 Necessity of extrema 

Multiplying equation by W* and integrating between 
[ri,r2] assuming, consistent with a sharp exponential de- 
cay, that W = or dW/dr — at these boundaries we find 
that 

r?l(£\\W\ a dr= rrnW\ 2 dr+ H ^\W'\ 2 dr 

(26) 

where derivatives are indicated with a prime. Since the right 
hand side is real, the imaginary part of the left and side must 
vanish. For general complex a this implies that 

7 I" < Z J -)'\W\ 2 dr = 0. (27) 

J ri {or + mil) 2 + Y \ 77 J 

Thus for a growing mode (7 7^ 0) to exist we need (c 2 /r/)' = 
at some r in [r%, r2\. For our disc, c 2 varies on a scale O(r), 
but r\ varies on a scale H <g r, thus given that the range 
of relative variation of the vortensity r) is of order unity, we 
infer that this quantity needs to have stationary points in 
order for there to be uns table modes. In the barotropic case 
IPapaloizou fc Linl l| 19891 ) have shown that vortensity max- 
ima are stable while minima are associated with instabili- 
ties. We expect these conclusions to apply here also because 
of the local nature of the modes of interest. We can ap- 
proximate c 2 ~ constant, or equivalently adopt a barotropic 
equation of state locally without changing the character of 
the problem. Referring back to Fig. |8(a)| it is clear that our 
basic state satisfies the necessary criterion for instability. 



4.3 Numerical solution of the eigenvalue problem 

We have solved the eigenvalue problem for the full equa- 
tion (|24[) using a shooting method that employs an adap- 
tive Runge- Kutta integrator a nd a multi-dimensional New- 
ton method (|Press et al.lll992T l. For low m 3) , unstable 
modes mainly comprise an evanescent disturbance near co- 
rotation (or the vortensity minimum at r = ro) and the 
simple boundary condition W = applied at the inner 
boundary r/r p = 0.55 and the outer boundary r/r p — 1.5 
produces good results. As m increases, the Lindblad reso- 
nances approach r — ro and a significant portion of the mode 
is wave-like requiring the application of outgoing radiation 
boundary conditions. We determined these using the WK BJ 
approximation (see eg. iKorvcanskv fc Papaloizoul (|l995l )V 

4-3.1 Eigenmode calculations 

We now discuss some example solutions to illustrate the in- 
stability of gap edges. We recall that the simulations indicate 
the ultimate dominance of small m values. One class of mode 
is associated with the inner vortensity ring while another is 
associated with the outer ring. 

As a typical example of the behaviour that is found for 
low m, each type of eigenmode for h = 0.05 and m = 3 is 
shown in Fig. [9] The background vortensity profile is also 
shown. The instabilities (7 < 0) are associated with the 
vortensity minima at t he inner a nd ou ter gap edge, as has 
also been observed by iLi et al l (|2005l ) in simulations, the 
modes are evanescent around co-rotation and the vorten- 
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sity peaks behave like wall through which the instability 
scarcely penetrates. The mode decays away from r — ro. 
For m = 3 Lindblad resonances occur at rz,/ro = 0.76, 1.21 
from which waves travelling away from co-rotation are emit- 
ted. However, the oscillatory amplitude is at most ~ 20% 
of that at r = ro. Hence for low-m the dominant effect of 
the instability will be due to perturbations near co-rotation. 
Increasing m brings tl even closer to ro, waves then prop- 
agate through the gap. The growth time-scale of the inner 
mode with m = 3 is ~ 14Po,- The outer mode has a growth 
rate that is about three times faster. Very similar results 
and growth rates were found for m = 1 and m — 2. Since 
the instability grows on dynamical time-scales, we expect 
non-linear interaction of vortices to occur within few tens of 
orbits and to affect planet migration if the latter were also 
on similar or longer time-scales. 

After the onset of linear instability and the formation 
of several vortices, it has been observed that non-linear ef- 
fects cause them to even tually merge into a single vortex 
|de Val-Borro et al.1 120071 ). This eventually interacts with 
the planet. In the fiducial simulation, rapid migration begins 
at 55Po which is compatible with the characteristic growth 
times found from linear theory. 

Fig. [5] indicates that the outer edge is more unstable 
than inner edge. The vortensity peaks are of similar height, 
but the inner minimum in r/ is less pronounced than outer 
because of the background profile. In this sense the outer 
edge profile is more extreme, and hence less stable. 

For illustration, we show in Fig. [10] eigenfunctions for 
the m = 7 mode of the outer ring. The equivalent mode was 
not found for the inner edge because high-m are quenched. 
Radiative boundary conditions were adopted in this case. 
Although the WKBJ condition is the appropriate physical 
boundary condition , its application here is uncertain be- 
cause the boundaries cannot be considered 'far' from the 
gap. However, solutions are act ually not sensitive to bound- 
ary conditions as reported by Ide Val-Borro et al.l (|2007h . 
Note the two spikes in Su r and 8u$ at r/ro — 0.90, 1.09 
which correspond to Lindblad resonances. These are not sin- 
gularities as can be seen from W(r) which is smooth there; 
other eigenfunctions were calculated from the numerical so- 
lution for W and thus may be subject to numerical errors. 
We see that increasing m increases the amplitude in the 
wave- like regions of the mode, but the growth rate is smaller 
than for m = 3. As the instability operates on dynamical 
time-scales, low-m modes will dominate over high-m modes, 
particularly through non-linear evolution and interaction of 
the former. 

We have also examined solutions with h = 0.03, 0.04 
and 0.06. In general, as h is lowered, 7 becomes more neg- 
ative. As the disc temperature is lowered with h, there are 
stronger shocks, and larger modifications to the disc profile, 
or steeper gradients. Hence we expect edges to become more 
unstable. In the case of h = 0.03, the rings become unstable 
(vortices form) before they reach an approximately axisym- 
metric state. 

Having understood the origin of vortices, we proceed 
to study their effect on planet migration, in the type III 
regime. 




Figure 10. m = 7 eigenmodes associated with the outer edge, 
obtained with WKBJ boundary condition. The perturbed quan- 
tities are plotted as \W\ (solid, black), \5u r \ (dotted, blue), 
(dashed, green) and scaled by their maximum values in 
r = [1.1,3.0]. The background vortensity profile is also shown 
(dash-dot, red), ro = 2.3 is the radius of the outer vortensity 



minima. 



5 SIMULATIONS OF FAST MIGRATION 
DRIVEN BY VORTEX-PLANET 
INTERACTION 

We now present hydrodynamic simulations of disc-planet 
interactions where the planet is free to migrate, so that r p = 
r p (t). Vortices that form in the co-orbital region near the gap 
edge move through the co-orbital region and cause torques 
to be exerted on the planet. The interaction between such 
large-scale structures and the planet leads to non-monotonic 
type III migration and is the focus of this section. The type 
III torque increases with the co-orbital mass deficit 8m: 



5m = 8-KTpBp 



s w(—x s ) 



w(x)dx 



(28) 



= S/w and Bp — 



( Mas set fc Papaloizou! I2003T I. where w 
^d r (r 2 Q) is the Oort constant evaluated at the planet ra- 
dius r v . Hence the inverse vortensity E/u> will be central to 
the discussion. 

The simulations described below have computational 
domain r = [ 0.4, 4.0 ] with resolution N r x N4, = 768 x 2304. 
The planet is set in circular orbit after which the initial ra- 
dial velocity is set t o be v r = —2>v/2r. as expec ted for a 
steady accretion disc |Lvnden-Bell fc P ringlc 1974|) The ini- 
tial surface density is chosen to be Eo = 7.0, corresponding 
to a few times the value appropriate to the minimum mass 
solar nebula in order to achieve rapid migration when a typ- 
ical viscosity vp = 1 is used |Masset fc Papaloizou! 120031 ; 
Ide Val-Borro et alj|2007h . For most of these simulations the 
full planet potential is applied from t = 0. Similar results 
were obtained if the potential is switched on over 5 orbits, 
which was employed in Section [3] when discussing formation 
of vortensity rings. In those cases, vortices were observed to 
form. Switching on the planet potential over several orbits 
does not weaken the instability, vortex-planet interactions 
still occur (see below), although at a slightly later time than 
if the planet is introduced at t — 0. 

Type III migration is numerical challenging due to its 
dependence on flow near the planet, one is sue being the 
numerical resolution. iD'Angelo et al.l (|2005T ) reported the 
suppression of type III migration in high resolution sim- 
ulations. The main migration feature discussed below is 
brief phases of rapid migration due to vortex-planet inter- 
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action, which does not depend on conditions very close to 
the planet. Our preliminary experiments with resolutions of 
N r x N<j, = 192 x 576, 256 x 768 both show such behaviour, 
thus we believe the higher resolution used below is sufficient 
to study this interaction. 

The locally isothermal equation of state implies no tem- 
perature changes due to the planet. This could lead to mass 
accumulation in the Hill sph e re and thus spurious torques 
from within. iPepliriski et al.l (|2008al ) used an equation of 
state where temperature increases close to the planet. We 
have also considered this model where the modification is 
applied to within the Hill sphere and at t = 0, c s (r = r p ) is 
18% higher than the local isothermal value. This again yields 
vortex-planet scattering, which is an indication that signif- 
icant torque contribution during the interaction is not due 
to material inside the Hill sphere. Qualitatively, the same 
phenomenon is observed for the case where the temperature 
increase is 38%. 

Finally, there is the issue of softening. We set e — 0.6H 
as before, but have considered softenings of 0.5H and 0.7H 
and both cases display vortex-planet interaction, although at 
earlier times when e is lower. This is expected since smaller 
softening has the same effect as increasing the planet mass, 
and thus producing a more unstable disc. 



5.1 Viscosity and vortices 

We have shown above that the ring structures formed by a 
Saturn mass planet in our case are linearly unstable. It has 
been shown through numerical simulations that su ch insta- 
biliti es are expected to lead to vortex production (|Li et alj 
l200ll ) and this has been seen when steep surface density 
gradients are pres ent outside the co-orbital regio n of a 



lower mass planet (jKoller et alj|2003l ; iLi et al 



monotonic migration was observed bv lOu et al 



2005). Non- 



(2007) when 



a vortex is present inside the co-orbital region, but the role 
of the vortex was not analysed in detail. 

Fixed-orbit simulations of disc-planet interactions show 
that a dimensi onless viscosity of order 1 0~ 5 suppresses vor- 
tex formation l|de Val-Borro et alj|2007l ). As a consequence 
studies using viscous discs typically yield smooth migration 
curves. We have confirmed this in our case with numerical 
simulations. 



5.2 Dependence of the migration rate on viscosity 

We first study type III migration as a function of viscos- 
ity. The effect of vortices appears at low viscosities. Fig. 1111 
shows the orbital semi-major axis a(t) for viscosities va = 
to vq = 1. As the orbit is very nearly circular, a(i) is always 
close to the instantaneous orbital radius r p (t). We consid- 
ered a case with vq — 10~ 3 which showed almost identical 
a(t) curve to vo = 0. The numerical viscosity is thus of or- 
der v = O(10 -8 ) which is much smaller than the typically 
adopted physical viscosity of v = 10~ 5 in such simulations. 

With the standard viscosity vo = 1, a — » a/2 in less 
than lOOPo, implying type III migration |Papaloizou et all 
120071 ). Comparing different vo, a(t) is indistinguishable for 
< t < 15, since viscous time-scales are much longer than 
the orbital time-scale. At t = 20, Idl increases with v§. In 



the limit v — > 0, the horse-shoe drag for a fixed orbit is oc v 
|Balmforth et al.ll200ll ; lMassetll2002j ). However, in this case 
d/0 there is a much larger rate-dependent torque respon- 
sible for type III migration, for which explicit dependence 
on viscosity has not been demonstrated analytically. 

Migration initially accelerates inwards (a, a < 0) and 
subsequently slows down at r ~ 1.4 (independent of vo). 
For vo = 1.0, 0.5 migration proceeds smoothly, decelerating 
towards the end of the simulation at which point the orbital 
radius has decreased by a factor of ~ 2.7 Migration curves 
for uq — 1.0 and vq = 0.5 are quantitatively similar. Low- 
ering i/o further enhances the deceleration at r ~ 1.4 until 
in the inviscid limit the migration stalls before eventually 
restarting. 

Despite differences in detail, the overall extent of the 
orbital decay in all of these cases is similar. This is expected 
in the model of type III migration where the torque is due to 
circulating fluid material switching from r p — x 3 — > r v + x s . 
In this model the extent of the orbital decay should not 
depend on the nature of the flow across r p , but only on the 
amount of disc material participating in the interaction, or 
equivalently the disc mass and this does not depend on v. 
On the other hand, the flow may not be a smooth function 
of time with migration proceeding through a series of fast 
and slow episodes as observed in Fig. 1111 Our argument is 
only valid if migration proceeds via the type III mechanism. 
We have not explored higher viscosities (y of more than a 
few times 10 ) in detail but migration time-scales of test 
cases indicate that in this regime, type III migration is not 
operating. 



5.3 Stalling of type III migration 

The issue discussed here is what inhibits the growth of |d|? 
Descriptions of (inward) type III migration usually assume 
that the libration time at r p —x s is muc h less than the time to 
migr ate across the co-orbital region |Masset fc Papaloizoul 
120031) . This implies that 



« 1, 



(29) 



where A p = 1/2(80./ dr) at r v (t) an d a is the changing semi- 
major axis. lPapaloizou et al.l l|2007h present a similar critical 
rate, but with the same dependence on Q, a, x s . If equa- 
tion (|29p holds, co-orbital material is trapped in libration 
on horse-shoe orbits and migrates with the planet. When 
X ^ 1 the horseshoe region shrinks to a tadpole, and mate- 
rial is trapped in libr ation about the L4 an d L5 Lagrange 
points (as observed bv lPeplihski et al"1l2008bf ). This can tend 
to remove the co-orbital mass deficit 8m which reduces the 
migration torque. Comparing \ f° r cases shown in Fig. 1121 
it is clear that migration with \ <C 1 ( hereafter in this 
context termed slow migration even though it may be much 
faster than type II migration) does not always hold, with 
max(x) ~ 0.6 being comparable for different vq. By follow- 
ing the evolution of a passive scalar we checked that horse- 
shoe material no longer migrates with the planet when |d| is 



1 We can regard the process of changing from a partial gap ex- 
tending for nearly the whole azimuth to one with a smaller az- 
imuthal extent, as gap filling. 
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during the planet modifies the co-orbital structure on orbital 
time-scales. Vortensity rings develop at r — i 
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Figure 11. Type III migration as a function of viscosity 




Figure 12. Evolution of Xi the ratio of libration-to-migration 
time-scale as a function of viscosity. Libration time is measured 
at r p — r = x s and migration time-scale is that across x s . 



large. This occurs for all v but only the low viscosity cases 
exhibit stalling. Hence, while horse-shoe material is lost due 
to fast migration, it is not responsible for stopping it0 By 
examining the inviscid case in detail later, we show that the 
stopping of migration is due to the flow of a vortex across 
the co-orbital region, where some of it becomes trapped in 
libration. 

5.4 The connection between the vortensity and 
fast migration 

The difference between the value of the inverse of the vorten- 
sity E/cj evaluated in the co-orbital region and the value 
associated with material that passes from one side of the co- 
orbita l region to the other d e fines the co-orbital mass deficit 
5m in Ma sset fc Papaloiz ou (2003). ft is often assumed that 
the vorticity is slowly varying so that the difference in the 
values of inverse vortensity reduces, to within a scaling fac- 
tor, simply to the difference in the values of the surface den- 
sity. 

Although Mas set &: Papaloizou! (|2003T ) assumed steady, 
slow migration in the low viscosity limit, it is nevertheless 
useful to examine its evolution in relation to the migration 
of the planet (Fig. III|) . Fig. [T3] shows the azimuthally av- 
eraged E/w-perturbation following planet migration. Intro- 



One may intuitively expect migration to be self-limited. 



~ ±2r h for 

t ^ 10 (Fig. 13(a) I. We showed above that this initial modi- 



fication is due to spiral shocks extending into the horse-shoe 
region and also that the ring-structure is unstable to the 
production of vortices. 

Increasing v reduces the rings' amplitude, but their 
locations are unaffected. Taking the length-scale of inter- 
est as I = rh ~ 0.1, for vo = 1 the viscous time-scale is 
td = I 2 /v — 56Po- hence at t < 10 viscous diffusion is not sig- 
nificant even locally. Thus ring-formation is not sensitive to 
the value of v. We note the correspondence between the sim- 
ilarity of the E/w profiles and similarity in a(t) for different 
v in the initial phase. That is, the co-orbital disc structure 
determines migration (|Masset fc Papaloizou! [20031 ) . Depen- 
dence on the value of the viscosity is seen beyond t = 50 
(Fig. 13(b) I, producing much smoother (and similar) pro- 
files for i/o = 0.5 and v = 1.0. 

For a fixed orbit, we would expect profiles to be 
smoothed on a viscous timescale. When there is migration, 
we should also consider advection of vortensity (in the planet 
frame). Following the planet, the perturbed profile can be 
smooth if the planet migrates through the background with- 
out carrying its original co-orbital material. This can be 
regarded as the loss of horseshoe material due to fast mi- 
gration. Any deviation from the background must then be 
due to local vortensity generation/destruction as the planet 
moves through (e.g. shocks). 

In Fig. 13(b) only the inviscid case is still in slow mi- 
gration, and only this disc retains the inner ring (low E/w). 
This suggests that the inner vortensity ring inhibits inward 
migration. In terms of 5m, for v 7^ the planet resides in a 
gap (co-orbital E/w is less than that at the inner separatrix, 
or 5m > 0) whereas in the inviscid case 5m ~ 0. 

Consider the v = case. The outer vortensity ring has 
widened to ~ 2r^ (c.f. Fig. 13(a) |. It is centred at 3r^ so 



that co-orbital dynamics may not account for it. However, 
the migration implies a flow of material across r p from the 
interior region. The increased region of low E/w exterior to 
the planet may be due to this flow. Notice the high-E/o; 
ring at +4r^ in Fig. 13(a) is no longer present in Fig. 13(b) 



because this ring is not co-orbital and therefore does not 
migrate with the planet. 



At t = 65 (Fig. 13(c) \ the vq — 0.5 and v = 1.0 cases 



continue smooth rapid migration (Fig. \\\\ with qualitatively 
unchanged profiles. For vq = 0.25, characteristic vortensity 
double-rings re-develop after a stalling event at t ~ 60 (Fig. 
1 1 1 P . The peaks and troughs of E/cj recover forms that are 



close to those in the initial phase (Fig. 13(a) |. At this time 
the inviscid case is in rapid migration. 

Fig. 1 1 3(d) I shows the final E/tj-perturbation profiles 
. Viscous cases are in slow migration, and have much 
smoother profiles. This can be due to diffusion (since we 
are at late stage of evolution) and/or migration across the 
background, in both situations there is little disc material 
carried by the planet. The inviscid case is also in slow mi- 
gration but retains the double-ring structure. This indicates 
two possible co-orbital configurations which slow down type 
III migration. 

In this section we reviewed the disc vortensity evolu- 
tion and its connection to the state of migration. A correla- 
tion between migration (fast, slow) and co-orbital structure 
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Figure 13. S/w-perturbation for different viscosities for the case 
So = 7 and M p = 2.8 X 10 -4 . Here we use u) = r~ 1 d r (ru^ > ) to 
approximate the vorticity. This is valid since \u r \ <C \uj,\. 



is apparent. The effective action of viscosity appears to be 
through its modification of the disc structure, rather than 
associated viscous torques acting on the co-orbital region. 



5.5 Evolution of the co-orbital region 

In Fig. [14] we illustrate the evolution of the following quan- 
tities associated with the co-orbital region that are related 
to the migration torque induced on the planet: 

(i) The gravitational torque acting on the planet due to 
fluid with |r — r p \ ^ 2.5rh, excluding fluid within of the 
planet. This includes co-orbital material and orbit-crossing 
fluid, the latter being responsible for the type III torque. 

(ii) The co-orbital mass deficit 8m is computed from az- 
imuthally averaged, one-dimensional disc profiles. We took 
the separatrices to be at \r — r p \ = 2.5rj,. 

(iii) Mtr, the mass of a passive scalar initially placed such 
that [r — r p \ = 2rn El- Note that Af tr = const, if it migrates 
with the planet. 

(iv) The average density (E) and vortensity (<^/E) of the 
region r — r p \ < 2.5rh- 

The evolution of viscous and inviscid cases is qualita- 
tively similar up to the stalling indicated by vertical lines in 
Fig. 1141 Prior to this there is a rapid migration phase asso- 
ciated with large a negative torque which we find originates 
from material crossing the planet orbit. Co-orbital torques 
are oscillatory and the period/amplitude is longer/larger for 
i^o = than for vo = 0.5. During rapid migration phases, the 
inviscid torque is twice as negative than in the viscous case. 
While the torque in a viscous disc remains negative, torques 
in an inviscid disc can be positive due to the formation of 



large-scale vortices. Note that the torque does not originate 
from within the Hill sphere since it is excluded from the 
summation. 

Migration is slow until sufficient difference builds up 
between co-orbital and circulating flow at which point there 
is a sudden flow-through the co-orbital region. At the same 
time there is significant loss of the original horse-shoe mate- 
rial (Af tr decreases by ~ 80%). The flow-through is reflected 
in vq = 0.5 (uo = 0) by a ~ 17% (36%) increase in (E) from 
t — 25 — > 50 (t = 50 — > 75). In the case of zero viscosity, this 
material is in the form of a high surface density vortex. We 



note from Fig. 14(b) that there is a repeated episode where 



there is a fall followed by a rise in (E) . However, in the late 
stages of evolution of the case with vq = 0.5, (E) decreases 
monotonically and the migration is slower. 

The co-orbital mass deficit is initially negativ^3 as the 
rings develop. It subsequently increases resulting in the onset 
of type III migration and is most positive during the follow- 
ing rapid migration phase, with peak values 8m x 10 3 ~ 1, 2 
for vq — 0.5, respectively. 8m x 10 3 then falls to +0.5 in 
the viscous case but to < in the inviscid case. In the lat- 
ter case, material flowing into the co-orbital region removes 
the co-orbital mass deficit and type III is suppressed; migra- 
tion experiences a more abrupt stall. 8m increases again for 
vo = while it remains approximately constant for Uq = 0.5 
and decreases towards the end. Type III migration can re- 
start in the inviscid disc but such behaviour is not observed 
for large viscosity. Type III is not operating in the late 
stages of the viscous case, in contrast to the inviscid evo- 
lution where fast type III migration is recurrent, faster than 
in the cases with applied viscosity, and is associated with 
large values of 8m. 

The above discussion shows that the magnitude of the 
applied viscosity is significant in determining the character 
of the migration. This is because the form of the flow through 
the co-orbital region is sensitive to the choice of viscosity. In 
particular, for high viscosity, this flow is smooth and there 
is less disruption of the co-orbital region. 



6 VORTEX-PLANET INTERACTION 

We now focus on the inviscid case where the role of vortices 
significantly affects the migration. We consider three phases 
apparent from (Fig. II 1 p . These are: a) t < 45 (slow migra- 
tion); &)60 < t < 70 (rapid migration); c) t ~ 75 (stalling) 
and d) 75 < t < 110 (second phase of slow migration). Typ- 
ical migration rates at various times are o(25) ~ —3 x 10 -3 ; 



d(65) 



-2 x 10~^ and d(85) ~-2x 10~ 3 . The vortex- 



induced rapid migration (phase b)) is almost an order of 
magnitude faster than the phases a) and d). Hence we refer 
to the latter as slow migration, but more accurately they 
are migration associated with gap formation. They are not 
necessarily slow in comparison to type I or type II migration. 

Different migration phases correspond to different disc 
structures. Fig. \T5\ shows global the global surface density 
evolution. At t — 24.75, during the slow migration phase 



3 Due to the uncertainty in the horse-shoe half-width x s we put 
the tracer well within the region defined by x s = 2.57"^ to ensure 
it is co-orbital. 



Ring structure in the vicinity of the separatrix means that the 
sign of <5m is sensitive to the adopted value of x s . Here we regard 
8m as a representation of gap depth, so we fix x s = 2.5r^. 
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Figure 14. Effect of viscosity on co-orbital evolution. The time 
evolution is illustrated for the model with v = 0.5 (upper panel) 
and for the inviscid model (lower panel). From top to bottom of 
each panel the evolution of the orbital radius, a(t), the co-orbital 
mass deficit 5m, the torque, the tracer mass Mti, the surface 
density (E) and the vortensity (^) are plotted. Angle brackets 
denote space averaging over the annulus [r p — 2.5ri l ,T p + 2.5r^ 
and A denotes perturbation relative to t = 0. The vertical line 
indicates stalling of migration. 
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Figure 15. The evolution of InE from early slow migration to 
stalling (~ 75 orbits). 



a), the planet resides in a partial gap (r ~ r p ± 2.3rh) 
with surface density ~ 20% lower than Eo. The gap is 
circular but non-axisymmetric with a low surface density 
arc along the outer gap edge trailing the planet. A par- 
tial gap is necessary for th e Type III migration mode 
l|Papaloizou fc Terqu cm 2006) but not sufficient. A surface 
density asymmetry ahead/behind the planet is neede d to 
provide the net co-orbital torque (| Artvmowic j l2004al ) . At 
this stage this is too weak. Further azimuthal density asym- 
metry has developed in the gap by t = 55 and there is a 
factor of ~ 3 variation in the gap surface density, but the 
asymmetry is still limited. 

Strong asymmetry can be provided by large-scale vor- 
tices near the gap edges of azimuthal extent A(j> ~ n (Fig. 
I15[) . The outer (inner) vortex rotates clockwise (counter- 
clockwise) relative to the planet so they are not co-orbital. 
Their origin was explained as a natural consequence of the 
instability of the vortensity rings (see section[4]above), their 
occurrence near gap edges has a strong influence on the type 
III migration mode. 

At t = 65.05 Fig. [T5] shows that the planet is just in- 
side the inner gap edge. The surface density contrast in the 
neighbourhood of the planet is largest as the inner vortex 
enters the co-orbital region from behind the planet. It ex- 
erts a net negative torque as the material crosses the planet 
orbit and enters (is scattered into) the exterior disc, and 
this snapshot corresponds to fast migration (phase b). At 
t = 75 the migration stalls and the planet no longer re- 
sides inside a gap. The planet has effectively left its gap by 
scattering vortex material outwards. This completes a single 
vortex-planet episode, during which the outer vortex simply 
circulates around the original outer gap edge and does not 
influence the co-orbital dynamics, though it contributes an 
oscillatory torque on the planet. 

The vortex-planet interaction is magnified in Fig. 1161 
At t = 65 the vortex circulates at ~ r p — 3r^ and has ra- 
dial extent ~ 3r^. The gap depth is largest and the mi- 
gration is fast. As the planet migrates inwards, the vor- 



Type III migration 15 







65.41 orbits, r„=1.62 /0.00 or;its, r p =1 A'i 




-6 -4 -2 2 4 6 
|r-r,!A 



Figure 16. Illustration of the surface density evolution from the 
start of rapid migration to stalling at t = 75. Maps of InS are 
plotted. 



Figure 17. Vortcnsity evolution: maps of the inverse vortensity 
ln(S/aj) are plotted spanning the time interval from early slow 
migration to stalling (~ 75 orbits) to the second phase of slow 
migration. 



tex splits with some material entering the co-orbital region 
while the rest continues to circulate (t = 70). Vortex ma- 
terial becomes trapped just behind the planet at t — 75, 
which would suggest a negative torque. However, the sur- 
face density distribution does not support the usual type III 
migration where horseshoe material moves with planet. Fur- 
thermore, the planet no longer resides in a gap. Part of the 
original horseshoe material is replaced with vortex material 
and gap-filling takes place. The co-orbital mass deficit is lost 
and the migration stalls. 

It is useful to examine the evolution of E/w, or inverse 
vortensity, since it define s co-orbital mass deficit 5m that 
drives the type III torque (|Masset fc Pa paloizou 2003). This 
is illustrated in Fig. 1171 In inviscid discs E/cj is approxi- 
mately conserved following a fluid so we can track material. 
The vortensity ring basic state and its stability was dis- 
cussed in section U (see Fig. [7}. By t — 55 the inner vortex 
has formed via non-linear evolution of the instability and 
begins to interact with co-orbital region. Vortensity conser- 
vation implies that the 'red blobs' near the outer ring were 
part of the inner vortex, consistent with inward migration of 
the planet via Type III (see Fig. I17[l . Each time this vortex 
passes by the planet, some of its material crosses the planet 
orbit from behind, thereby exerting a negative torque. Rapid 
migration at t = 65 occurs when the main vortex body flows 
across the co-orbital region. The inner ring is disrupted and 
no longer extends 2tt in azimuth. 

This contrasts with the usual type III scenario where 
material simply transfers from inner to outer disc leaving 
the co-orbital region unaffected. In the inviscid disc, dis- 
ruption is necessary due to the existence of vortensity rings 
of much higher vortensity than the vortex. Under type III 
and vortensity conservation, vortex material must cross the 
planet orbit without changing its vortensity. This would not 
be possible if a ring structure is maintained. In this sense, 
the vortensity rings oppose the type III mode. Hence, mi- 



gration is slow until significant ring disruption occurs that 
is associated with the vortex flowing across. 

When migration stalls at t = 75, Fig. [17] shows that 
the vortensity rings are much less pronounced compared to 
initial phase. The vortex splits into several smaller patches 
circulating in the original gap. At the planet's new radius, 
material of high E/oj fills the new co-orbital region, corre- 
sponding to lower co-orbital mass deficit. However, by t — 85 
(not shown) new vortensity rings are setup near the new or- 
bital radius and are qualitatively similar to those present in 
the ring basic state. The vortex material which passed to 
the outer disc during the first rapid migration phase is now 
irrelevant, much like the outermost vortex. During the build 
up to the second rapid migration phase at t = 100, there is 
a single vortex associated with inner gap edge and two asso- 
ciated with the outer ring. This simply leads to a repeat of 
the first rapid migration phase. However, the presence of a 
vortex inside the co-orbital region, from the first rapid mi- 



gration phase reduces the co-orbital mass deficit (Fig. 14(b) I 



and hence the migration rate for the second rapid phase (at 
* ~ 140). 

6.1 The effect of changing the disc mass 

Here we consider type III migration in inviscid discs of dif- 
ferent masses obtained by scaling the parameter (Eo). Sim- 
ulations here had a resolution N r x N^, — 512 x 1536. 

Fig. [18] shows migration as a function of Eo For Eo = 
5 — 9 the planet migrates by the same amount during the 
first rapid migration phase and stalls at the same radius. 
This is also observed for a second rapid migration phase for 
Eo = 6 — 9. Although the highest density case is unphysical 
due to the lack of inclusion of self-gravity, results are con- 
sistent with the notion that rapid migration is initiated by 
sufficient contrast between co-orbital (gap) and circulating 
fluid (vortex), measured by 5m. 

The jump in orbital radius, should it occur, is indepen- 
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Figure 18. Vortex-induced migration in discs with different scal- 
ings of the initial surface density profile. The extent of rapid mi- 
gration is independent of the initial surface density scaling. 



dent of So- As the interaction involves the vortex flowing 
from the gap edge across the co-orbital region; its change 
in specific angular momentum is independent of density, be- 
cause the co-orbital region size is fixed by planet mass. The 
results then suggest the vortex needs to grow to a mass M v , 
only dependent on the planet mass, in order to scatter the 
planet. Since the vortex forms at the gap edge, M v can be 
linked to Sm because the co-orbital mass deficit depends on 
the edge surface density. 

As the vortex originates from instabilities with growth 
rate independent of So, increasing So means less time is 
needed for the vortex to build up to critical mass or density. 
Hence, increasing surface density only shortens the 'waiting 
time' before rapid migration. However, if the density is too 
low, e.g. So = 3 then vortex-induced rapid migration may 
never occur. 

Consider an inner vortex of mass M v formed by insta- 
bility at r v — r p — fir^ (f3 > 0) with width ar^, where the 
Hill radius r h = f r p , fo = (M P /(3M*)) 1/3 . M v is clearly 
limited by the amount of material that can be gather into 
the vortex, so that M v < 2irarhr v Y,. Taking M v /M p = 3.5 
as critical for rapid migration 0, this means 



So > 



3.5M„ 



27ra/o(l-/3/ )r| 



x 10 4 



(30) 



Taking representative values found in simulations of a = 
3, f3 = 4 and r p ~ 2 gives So > 3.5. For such cases rapid, 
vortex-induced migration was indeed observed fFig |18[) but 
for So = 3 it was not. This is similar to the usual require- 
ment that in order for Type III migration to occur for in- 
termediate planets, the disc should be sufficiently massive 
l|Masset fc Papaloizou! [20031 ) . In our case the limitation is 
specifically due to the maximum possible vortex mass. 



6.1.1 Critical co- orbital mass deficit 

In order to link type III migration and vortex-planet scat- 
tering, we measured the co-orbital mass deficit Sm, which 
amounts to comparing average inverse vortensity of co- 
orbital fluid just behind the planet, to that of circulating 



5 We estimate the vortex mass by monitoring the decrease in 
total mass of an annulus interior to the inner gap edge. 




(c) So 



(d) So 



Figure 19. Evolution of the co-orbital mass deficit as defined 
by inverse vortensity. The average S/oj for co-orbital material 
is taken over r — r p = [— 2.5r^, 0], <j> = [<j> P ,(p P — Jr/4]; and that 
of circulating material is taken over r — r p = [— 6, — 2.5]r h , <j> = 
[4> p , <f> p — 7r /4] . Very similar behaviour was obtained when we com- 
pared the co-orbital and circulating surface densities. 



fluid just inside the inner separatrix but also behind the 
planet: 

Sm = 27rx a r~ 1/ ' 2 ((S/a;) c i rc — (S/tj) coor b). 

This is a simplified v ersion of the definition in 
iMasset fc Papaloizou! i|2003ft . Results are shown in Fig. 1191 
The oscillatory nature of 8m reflects a vortex circulating at 
the gap edge, Sm maximising when the vortex is within the 
patch of fluid where averaging is done. As it grows, the high 
vortensity vortex contributes to (S/cj) c i rc hence favouring 
type III migration. Cases with rapid migration share the 
same evolution of 8m. 8m increases up to ~ 4-5 before the 
vortex first induces fast migration. For So = 3, which does 
not show such a rapid migration phase, typically Sm < 5 
with smaller amplitude variation. 

During fast migration Sm rapidly decreases and migra- 
tion stalls when Sm < 0. This is because as the vortex flows 
across the co-orbital radius, it contributes to (S/oj) coor b, 
lowering Sm. However, we comment that details depend on 
the unperturbed (t — 0) surface density profile. Discs dis- 
cussed here initially have uniform surface density. If a sur- 
face density S oc r~ p , p > were adopted, the planet can 
be scattered to a region of higher background S/cj com- 
pared to the flat case, so (S/ui} c irc may increase due to the 
background. We have run simulations with p — 0.3, 0.5 and 
found that the periods of stalling are of shorter duration, 
consistent with the discussion above regarding the variation 
with the surface density scale. 



6.2 Comparison to Jupiter-mass planet 

For completeness we briefly consider the case of a Jupiter- 
mass planet. The setup is the same as the previous sub- 
section (with So = 7), but the planet potential is now 
switched on over 5 orbits since Jupiter is 3 times more mas- 
sive than Saturn. Fig. [20] compares the migration curves 
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a{t) for M p = 10~ 3 in inviscid and viscous (vo = 1, 2) 
discs to a Saturn mass planet in a inviscid disc. The same 
vortex-planet scattering occurs for the Jupiter-mass inviscid 
case. This was checked explicitly by examining the vorten- 
sity evolution and observing a vortex associated with the 
inner gap edge build-up and flowing across the co-orbital re- 
gion. Jupiter induces stronger shocks and therefore higher- 
amplitude vortensity rings, so they are more unstable. The 
instability growth time-scale is therefore shorter than that 
of Saturn, so the vortex-induced migration occurs very soon 
after the planet in introduced. 

The a(t) for Jupiter in an inviscid disc show two mi- 
gration phases — fast and slow, and that significant or- 
bital decay occ ur within ~ 50Pp. T hese features were 
also obs erved b y Pcplihs ki et al.l (|2008bT ). However, the fast 
phase of lPeplihski et all is almost a linear function of time, 
whereas in our inviscid case migration is clearly accel- 
erating inwards during t < 50Po. In our case there is 
an abrupt tran sition to the slow phase whereas that of 
IPeplihski et alj is smooth er and no vortices were identified. 
Although iPeplihski et ail did not include physical viscosity, 
their Cartesia n code is more diffusive a nd vortex-formation 
is suppressed |de Val-Borro et al.l [20071. In the late stages 
of our Jupiter in the inviscid disc, thin vortensity rings re- 
form and inhibit further migration. A new vortex develops 
at this stage but it does not grow a sufficient size to induce 
another episode within simulation time. 

Interestingly, increasing to standard viscosity i/ = 1 
results in a similar behaviour for the fast phase, though the 
transition to slow migration is smoother than zero viscos- 
ity. The vortensity distribution for vo = 1 is also much 
smoother than v§ = and individual vortices were not 
identified for the inner g ap edge. This is consistent with 
Ide Val-Borro et all (|2007h who showed that v = O(10 -5 ) 
suppress vortex formation. With Uq = 1, flow-through the 
co-orbital region is a smoother function of time, unlike 
episodic behaviour of inviscid discs, and therefore do not 
experience sudden stalling. 

A case with i/o — 2 is also show in Fig. 1201 The orbital 
decay time-scale is again ~ 50Po, consistent with type III 
migration, although there could be complications from the 
fact that a Jovian-mass planet in a viscous disc can undergo 
type II migration. The transition from fast to slow is again 
smoothed (cf. vq = 0, 1), the fast phase is no w more linear, 
qualitatively closer to IPeplihski et all |2008b). These results 
suggest that our invi scid cases have l ower total viscosity than 
those considered by IPeplihski et all 



7 SUMMARY AND DISCUSSION 

In this paper we have studied the role of large-scale vor- 
tices in type III migration in low viscosity discs. We focused 
mainly on Saturn-mass planets because they open partial 
gaps, a configuration where type III migration c an operate 
if the disc is massive (|Masset fc Papaloizoull2003l ). Type III 
migration would occur when the planet is Saturn mass be- 
fore growing to Jovian mass. We first demonstrated through 
numerical simulations and semi-analytic modelling of invis- 
cid discs, that vortensity rings originate from spiral shocks 
induced by the planet. For Saturn or more massive planets, 
the rings reside just inside the co-orbital region, while for less 




Figure 20. Vortex-induced migration for a Jupiter-mass planet 
in an inviscid disc (dotted) and viscous disc (vq = 1, dashed; 
vq = 2, dashed-dot) compared to Saturn-mass planet in a inviscid 
disc (solid). 



massive planets they are not co-orbital features. Vortensity 
rings are setup whether the planet is introduced suddenly or 
switched on over several orbits and instability occur in bo th 
approaches (the latter used by Ide Val-Borro et alj (|2007l )'). 

We showed that the gap edges, associated with lo- 
cal vortensity minima, are dynamically unstable to non- 
axisymmetric perturbations. Dominant unstable modes are 
localised go on to develop into vortices which merge in the 
non-linear regime, as verified by simulations. The case of a 
Jupiter-mas s planet held on a fixed o rbit has already been 
discussed by Ide Val-Borro et alj l|2007l ) which shows vortex- 
formation at gap edges. The effect of vortices on migration 
is most significant in low viscosity discs (y < 0.25 x 10 ) 
because the instability is suppressed at higher viscosity. In 
the inviscid limit only a small numerical viscosity is present. 
Suppose we consider a physical viscosity of v — O(10 -6 ) , 
it is still large compared to the numerical viscosity, so we 
expect the latter to be unimportant. However, v — O(10 -6 ) 
is already unable to suppress the instability. Our conclu- 
sions are unaffected by numerical viscosity of this code. The 
presence of high density vortices at the gap edge produces 
non-smooth migration, with episodes of fast migration corre- 
sponding to the vortex-planet interaction. This is analogous 
to planet-planet scattering, and the planet's orbital radius 
jumps by a few Hill radii in one episode. The vortex is also 
responsible for stalling migration in discs with initially flat 
surface density. In this case there can be repeated episodes of 
vortex-induced migration. Viscosity smooths the flow across 
the co-orbital radius, but has limited effect on the extent of 
orbital decay via the type III mechanism. 

We also explored the role of vortices in inviscid discs 
of different masses. The extent of orbital decay in a single 
episode of fast migration is independent of initial surface 
density scaling. This suggests a critical vortex mass or sur- 
face density is required to interact, which can be linked to 
the concept of the co-orbital mass deficit that drives type III 
migration. The case of a Jupiter-mass planet in an inviscid 
disc also displayed vortex-planet interaction. 



7.1 Outstanding issues 

One issue not considered in detail in this paper are bound- 
ary effects, although we used the damping boundary con- 
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ditions l|de Val- Borrol l2006h which is aimed to remove re- 
flections. The instabilities that lead to vortex development 
are localised near gap edges. As long as the planet, and 
hence vortex-formation, is far from the disc edges, bound- 
aries effects on vortex-planet interaction should be limited. 
We have run a simulation where the domain is expanded to 
r = [0.2, 6.0] and compared a(t) to the fiducial disc which 
extends r = [0.4,4.0]. The extended domain case still ex- 
hibited non-smooth a(t) with similar orbital decays during 
the first two episodes of rapid migration. As the domain was 
extended inwards, there was a third phase of rapid migra- 
tion. There is more discrepancy later in the evolution as the 
planet is migrating towards the inner boundary. Boundary 
effects need to be examined in more detail in the future, but 
should not change the main conclusion of this paper, which 
is that in a low viscosity disc, type III migration is induced 
by vortices and episodic. 

Although we have demonstrated the effect of vortices 
on migration, it is natural to question their existence as 
long-li ved struc t ures i n real discs with finite viscosity. Re- 
cently, iLi et all (|20Qgh studied migration in nearly-laminar 
discs in the type I regime. They found effects at low viscos- 
ity, which tend to slow down migration, begin to appear for 
v ^ 10~ 7 . Their simulations last O(10 3 ) orbits. In our case, 
non-smooth migration can be observed at v — O(10 -6 ). We 
adopted a more massive planet so vortex formation occurs 
quickly, and since type III migration time-scales are much 
shorter than type I, a very small viscosity is not required 
in order to allow the development of sharp features in the 
disc (which the instability relies on). Vortices can thus inter- 
fere with migration when viscosity is one order of magnitude 
smaller than the standard value v = 10 -5 . 

We have used the simplest model and numerical setup 
to describe the disc-planet system. One physical issue is 
the lack of inclusion of self-gravity (SG). It may be im- 
portant when discussin g type III migration since this op- 
erates in massive discs (|Masset fc Papaloizou! l2~003h . In or- 
der to ha ve self-consisten t physics, self-gravity is essential. 
Although ILi et all ((2009) included SG, its role was not dis- 
cussed. The effect of self-gravitating vortices on the migra- 
tion of the larger mass planets considered here should be ex- 
plored. We note in the standard viscous disc with v — 10~ 5 , 
vortices are transient features. However, it is conceivable 
that SG may mitigate the effects of viscous diffusion be- 
cause vortices have high surface density. Thus, vortices may 
exist for higher viscosity values in SG discs. These issues will 
be the subject of a future study. 
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APPENDIX A: UPPER LIMIT ON THE 
HORSESHOE WIDTH 



APPENDIX B: VORTICITY JUMP ACROSS A 
STEADY ISOTHERMAL SHOCK 



We deduce an upper limit on the horseshoe width x 3 that 
is valid in the limit of either zero or constant pressure. Con- 
sider a local Cartesian frame that has origin at the planet 
and hence co-rotates with it with the Keplerian angular ve- 
locity Q p — yj GM, I r-3 . Here M* is the mass of the central 
object and r p is the orbital radius of the planet. In this frame 
fluid elements approach the planet from (x — xq > 0, y = 
oo) with velocity (0, — ZQ.pXo/2). Suppose such a fluid ele- 
ment executes a horseshoe turn crossing the co-orbital radius 
at q = (0, y) with velocity (v x , 0). 

We wish to determine an upper bound on the value of 
xq for which such motion occurs. The equation of motion 
governing a fluid element or particle is 



Dv 



Dt 
where 



+ 2fi p zAv = -V$cff, 



GM V 



^Q. 2 x 2 



\J x 2 + y 2 2 p 



(Al) 



(A2) 



Here the effective potential $ c ff contains contributions from 
the gravitational potential of the planet of mass M p and the 
tidal potential associated with the central object. 

From equation (|A1[) it follows that the Jacobi invariant 



J = \f + $cft 



(A3) 



is constant along a particle path. Equating J evaluated at 
(x,y) — (xq,oo) to J evaluated at q gives 



o 2 2 — 2 

-ilpXQ — — V x 



V 



(A4) 



The steady-state Euler equation of motion for v y evaluated 
at q is 



v x d x Vy + 2Q. p v x = —GMp/y 2 



(A5) 



Since in the neighbourhood of q for the type of streamline 
we consider, v x < and d x v y < we deduce that 



I . GMp 



(A6) 



Combining equation (|A4[I and equation (|A6f) we find that 



8 p 0> 2 {2n p y 2 



GM V 



(A7) 



Writing xo = Xorh, where the Hill radius r% = 
(M p /(3M*)) 1,/3 rp and similarly setting y = yrn, we obtain 



xo < 



8 _ _3 

y V 4 



< 2.3. 



(A8) 



Thus because the maximum possible value of the right hand 
side of the above is 2.3, we deduce that particles execut- 
ing a U-turn could not have originated further than 2.3rh- 
This is comparable to the value of 2.5r^ that has been esti- 
mated from hydrodynamic simula tions (|Artvmowic z 2004b; 
iPaardekooper fc Papaloiz"oull2009l '). 



Consider an isothermal shock that is stationary in a frame 
rotating with angular velocity fl p z. In order to evaluate the 
vorticity jump across the shock it is convenient to use a 
two dimensional orthogonal coordinate system (3:1,2:2) de- 
fined in the disc mid-plane such that one of the curves 
X2 = constant = x s coincides with the shock. The curves 
x\ = constant will then be normal to the shock where they 
intersect it. In addition the coordinates are set up so that 
(xi,X2,z) is a right handed system. The z-component of 
relative vorticity uj r can then be written as 



h\hi 



d{u2h,2) d(uihi) 



d. 



8X2 



(Bl) 



where (fti,/i2) are the components of the coordinate scale 
factor. 

We note that on X2 — x s , ui is the velocity tangential 
to the shock. The normal component 112 and other state 
variables undergo a jump from pre-shock values to post- 
shock values on normally traversing the curve X2 = x s . For 
an isothermal shock 

^E^L =M - 2 = (B2) 

U2 Ppost 

where M = U^/d, is the pre-shock perpendicular Mach num- 
ber, and here and in similar expressions below connecting 
pre-shock and post-shock quantities, we have denoted post- 
shock values with a subscript post while leaving pre-shock 
quantities without a corresponding subscript. Thus the jump 
in relative vorticity is 



(B3) 



Quite generally the xi component of the equation of 
motion for a steady state flow is 



ui du\ U2 dui f U2 dli2 Mi dhi 

hi dxi h,2 dx2 \h1h2 dxi h\li2 8x2 

1 dP 19$ 

= T - q 1" 2S2pU 2 , 

phi 0x1 hi 0x1 

or equivalent ly 

iti dui u 2 du 2 , . 

TT> '"TTa (2ilp + uJr)u 2 = 

hi 0x1 hi 0x1 phi 0x1 



(B4) 



1 bp 1 a$ 

hi dxi ' 
(B5) 



where P is the pressure and $ the total potential (the latter 
quantity being continuous across the shock). From (|B5|) we 
obtain an expression for the relative vorticity in the form 



1 du2 ui dui 1 

uj r = 1 1 

hi dxi u-ihi dxi pU2hi dxi u^hi dxi 



dP +^^-2n v 



(B6) 

Applying equation (|B6|) to give an expression for the post 
shock relative vorticity, and using equation (|B2() to express 
the post-shock normal velocity and density in terms of the 
corresponding pre-shock quantities, the vorticity jump can 
be written in the form 

. 1 d{M- 2 u 2 ) M 2 m dui M 2 <9$ 
|Wrl — 1 ; — 7; V 



+ 



hi dxi U2hi dxi U2/11 dxi 

1 <9-Ppost 

pu2hi dxi 



(oj r + 2Q,p) 



(B7) 
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Adopting a locally isothermal equation of state, we have 
fpost = M 2 P. Substituting this into equation (|B7|) . while 



making use of equation (|B6[) together with the relation U2 
c s M, we obtain 



M 2 In dxi 
(M 4 - 1) dc s 



(B8) 



Mhi dxi ' 

- uj r is the absolute vorticity. 
Setting hidxi = dS with dS being the corresponding 
element of distance measured parallel to the shock, this takes 
the form 

r i r i C S (M 2 - l) 2 3M ,„, 2 

(M 4 - 1) 0c s 

M~ dS' (B9) 

This gives the vorticity jump across a shock in terms of pre- 
shock quantities measured in the rotating frame in which it 
appears steady. It is important to note that the expression 
()B9|) applies specifically in the right handed coordinate sys- 
tem we have adopted with shock location given by X2 — x s . 
If we instead adopt xi = x s for this location, the signs of 
the derivative terms wou l d be reversed as in the expression 
(2.23) given in lKevlahanl (|l997l ). Note too that the last term 
on the right hand side that is proportional to the gradient 
of c s along the shock arises from the assumption of a lo- 
cally isothermal eq uation of s t ate a nd is not present in the 
treatment given bv lKevlahanl(ll997l ). 



